Lecture 10: Cycles

1. Cyclostratigraphy
2. How to identify a cycle
3. Fourier Transforms

We acknowledge and respect the lək̓ʷəŋən peoples on whose traditional territory the university stands and the Songhees, Esquimalt and W̱SÁNEĆ peoples whose historical relationships with the land continue to this day.

Cyclostratigraphy

• Cyclostratigraphy is a sub-discipline of stratigraphy that seeks to identify, characterize and interpret cyclic variations in the stratigraphic record
  • identify cycles $\rightarrow$ interpret timing of cycles $\rightarrow$ age models and correlations
  • fundamentally, it is about explaining the processes behind a record

Identifying cycles

from scipy.fftpack import fft
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
%matplotlib inline

def plot_data_res(x,y):
    fig1=plt.figure(1)
    fig1.clf()
    ax1=fig1.add_subplot(121)
    
    ax1.plot(x,y,'-',color='#1E90FF')
    
    z=np.polyfit(x,y,1)
    p=np.poly1d(z)
    
    ax1.plot(x,p(x),'--',color='#B22222')
    ax1.text(0.05,0.91,r'r$^2$ = %2.2f' %(np.corrcoef(x,y)[0,1]**2),
             horizontalalignment='left',transform=ax1.transAxes)
    
    ax2=fig1.add_subplot(122)
    res=y-p(x)
    
    ax2.hist(res,color='#808080')
             
    ax2.text(0.05,0.91,r'$\mu$ = %2.1f' %(np.mean(res)),
             horizontalalignment='left',transform=ax2.transAxes)
    ax2.text(0.05,0.85,r'$\sigma$ = %2.1f' %(np.std(res)),
             horizontalalignment='left',transform=ax2.transAxes)

Identifying cycles

def plot_data_res_fft(x,y,yf,xf):
    fig1=plt.figure(1)
    fig1.clf()
    ax1=fig1.add_subplot(221)
    ax1.plot(x,y,'-',color='#1E90FF')
    z=np.polyfit(x,y,1)
    p=np.poly1d(z)
    ax1.plot(x,p(x),'--',color='#B22222')
    ax1.text(0.05,0.91,r'r$^2$ = %2.2f' %(np.corrcoef(x,y)[0,1]**2),
             horizontalalignment='left',transform=ax1.transAxes)
    ax2=fig1.add_subplot(222)
    res=y-p(x)
    ax2.hist(res,color='#808080')      
    ax2.text(0.05,0.91,r'$\mu$ = %2.1f' %(np.mean(res)),
             horizontalalignment='left',transform=ax2.transAxes)
    ax2.text(0.05,0.85,r'$\sigma$ = %2.1f' %(np.std(res)),
             horizontalalignment='left',transform=ax2.transAxes)
    ax3=fig1.add_subplot(212)
    ax3.plot(xf,2.0/N *np.abs(yf[0:N//2]))
    ax3.text(0.99,0.9,'a periodogram',transform=ax3.transAxes,horizontalalignment='right')
    
    power=2.0/N *np.abs(yf[0:N//2])
    
    max_power=np.max(2.0/N *np.abs(yf[0:N//2]))
    freq_idx=np.where(power==max_power)[0]
    freq=np.round(xf[freq_idx],2)
    
    ax3.plot(xf[freq_idx],max_power,'ro')
    
    return [ax1,ax2,ax3],max_power,freq  

Identifying cycles

t = 20 # total time series
N=200 # sample spacing 
x = np.linspace(0.0, t, N) #noise amplitude
a=2
noise=np.random.normal(0,a,N) #total signal
y=x+noise # plot_data_res(x,y)

Identifying cycles

#add a sine wave
y=y+8*np.sin(0.5*2*np.pi*x)

# plot_data_res(x,y)

Identifying cycles

#remove the linear rise
y=signal.detrend(y)
#define sampling interval for frequency axis
dt=np.round(np.diff(x),10)[0]
#perform Fast Fourier transform  
yf = fft(y)
xf = np.linspace(0.0, 1.0/(2.0*dt), N//2)
#plot results
# axes,max_power,freq=plot_data_res_fft(x,y,yf,xf)

Identifying cycles

#create and plot the recovered sinusoid
amp=max_power
ang_freq=freq * 2*np.pi
phase=0 * np.pi/180
cycle=amp*np.sin(ang_freq*x + phase)
#axes,max_power,freq=plot_data_res_fft(x,y,yf,xf)
#axes[0].plot(x,cycle,'k--')

Identifying cycles

##recalculate residuals
#axes=plot_data_res_fft(x,y-cycle,yf,xf)

quickly we cannot rely on our eyes ...

quickly we cannot rely on our eyes ...

quickly we cannot rely on our eyes ...

quickly we cannot rely on our eyes ...

(Draw cycles in sediments)

Information can be described as an infinite series of sin and cos waves of differing frequencies and magnitudes.

Euler's formula establishes the fundamental relationship between the trigonometric functions and the complex exponential function.

$$ e^{ix} = cos x + i~sin x $$"the most remarkable formula in mathematics" - Richard Feynman

np.exp(1j*3) 

(-0.9899924966004454+0.1411200080598672j)

x=np.linspace(0,2*np.pi,100)
y=np.ones_like(x)
y[x>np.pi]=6
C=y*np.exp(1j*x/2)
plt.plot(C.real,C.imag,'.')

[<matplotlib.lines.Line2D at 0x7ffa97a744f0>]

plt.figure(figsize=(16,8))
plt.subplot(1,2,1)
xt = np.linspace(0, 4, 1000)
y = np.sin(xt*np.pi)
plt.plot(xt, y)
plt.subplot(1,2,2)
freq=1/2
complex_y = y * np.exp(-2 * np.pi * 1j * freq   * xt) #using imaginary j  and exp to convert to polar
plt.plot(complex_y.real, complex_y.imag)

[<matplotlib.lines.Line2D at 0x7ffa95c7fb80>]

plt.figure(figsize=(25,5))
plt.subplot(1,2,1)
xt = np.linspace(0,4,1000)
y = (np.cos(xt*3*2*np.pi)+1)/2
plt.plot(xt,y,lw=2)
plt.gca().set_title('3 beats/second',y=1.1)
plt.gca().set_xlabel('Time (s)') 
plt.subplot(1,2,2, polar=True)
_=plt.gca().set_yticklabels([])  
plt.plot(xt,y)

[<matplotlib.lines.Line2D at 0x7ffa95b0d370>]

plt.figure(figsize=(25,5))
plt.subplot(1,2,1)
xt = np.linspace(0,4,1000)
y = (np.cos(xt*3*2*np.pi)+1)/2
plt.plot(xt,y,lw=2)
plt.gca().set_title('3 beats/second',y=1.1)
plt.gca().set_xlabel('Time (s)')
plt.subplot(1,2,2, polar=True)
plt.plot(xt*.5*2*np.pi,y,alpha=1,lw=2)
plt.gca().set_yticklabels([])
_=plt.gca().set_title('0.5 cycles/second',y=1.1)

Some set up to cycle through different frequencies:

%%capture
import time
import pylab as pl
from IPython import display
xt = np.linspace(0,4,1000)
y = (np.cos(xt*3*2*np.pi)+1)/2
y = (np.cos(xt*3*2*np.pi)+1)/4+(np.cos(xt*2*2*np.pi)+1)/4
plt.gca().set_title('3 beats/second',y=1.1)
plt.gca().set_xlabel('Time (s)')

Wrapping the signal (left) around the unit circle (right) at different frequencies:

for i in np.linspace(.1,6,100):
    fig=plt.figure(figsize=(15,5))
    plt.subplot(1,2,1)
    plt.plot(xt,y,lw=2)
    plt.subplot(1,2,2, polar=True)
    display.clear_output(wait=True)
    plt.plot(xt*i*2*np.pi,y,alpha=1,lw=2)
    plt.gca().set_yticklabels([])
    _=plt.gca().set_title(f'{i:.2f} cycles/second',y=1.1) 
    display.display(pl.gcf())
    time.sleep(.25)

plt.figure(figsize=(25,5))
plt.subplot(1,2,1)
xt = np.linspace(0,4,1000)
y = (np.cos(xt*3*2*np.pi)+1)/2
plt.plot(xt,y,lw=2)
plt.gca().set_title('3 beats/second',y=1.1)
plt.gca().set_xlabel('Time (s)')
plt.subplot(1,2,2, polar=True)
# plt.plot(xt*.5*2*np.pi,y,alpha=1) 
plt.plot(xt*5*2*np.pi,y)
plt.gca().set_yticklabels([])
plt.gca().set_title('5 cycles/second',y=1.1)

Text(0.5, 1.1, '5 cycles/second')

plt.figure(figsize=(25,5))
plt.subplot(1,2,1)
xt = np.linspace(0,4,1000)
y = (np.cos(xt*3*2*np.pi)+1)/2
plt.plot(xt,y,lw=2)
plt.gca().set_title('3 beats/second',y=1.1)
plt.gca().set_xlabel('Time (s)')
plt.subplot(1,2,2, polar=True)
# plt.plot(xt*.5*2*np.pi,y,alpha=1) #add back
# plt.plot(xt*5*2*np.pi,y,alpha=1)
plt.plot(xt*3*2*np.pi,y)
plt.gca().set_yticklabels([])
plt.gca().set_title('3 cycles/second',y=1.1)

Text(0.5, 1.1, '3 cycles/second')

plt.figure(figsize=(25,5))
plt.subplot(1,2,1)
xt = np.linspace(0,4,1000)
y = (np.cos(xt * 3 * 2 * np.pi) + 1) / 4 + (np.cos(xt * 2 * 2 * np.pi) + 1) / 4 - .5
plt.plot(xt,y)
plt.gca().set_title('3 beats/second + 2 beats/second',y=1.1)
plt.gca().set_xlabel('Time (s)')
plt.subplot(1,2,2, polar=True)
plt.plot(2 * np.pi * 0.5 * xt,y,label='0.5 c/s')
plt.plot(2 * np.pi * 2 * xt,y,label='2 c/s')
plt.plot(2 * np.pi * 3 * xt,y,label='3 c/s')
plt.gca().set_yticklabels([])
plt.legend(loc='lower right',fontsize=10,bbox_to_anchor=(1.3, -.1))
# plt.gca().set_title('3 cycles/second',y=1.1)

<matplotlib.legend.Legend at 0x7f780e399fd0>

plt.figure(figsize=(25, 5))
plt.subplot(1, 2, 1)
xt = np.linspace(0, 4, 1000)
y = (np.cos(xt * 3 * 2 * np.pi) + 1) / 4 + (np.cos(xt * 2 * 2 * np.pi) + 1) / 4 - .5
complex_y = y * np.exp(-2 * np.pi * 1j * 3 * xt) #using imaginary j  and exp to convert to polar
plt.plot(complex_y.real, complex_y.imag)
plt.gca().set_aspect(1); plt.gca().set_title("3 cycles/second"); plt.subplot(1, 2, 2)
x_COM = []
freqs = np.linspace(0.1, 10, 1000)
# freqs = np.arange(0,1000-1)
for i in freqs:
    complex_y = y * np.exp(-2 * np.pi * 1j * i * xt)
    x_COM.append((np.sum(complex_y.real)**2+np.sum(complex_y.imag)**2)**(1/2))
     
plt.plot(freqs, x_COM,alpha=1,zorder=2); _ = plt.gca().set_xticks(range(0, 11)); plt.gca().set_ylabel("'sum'"); plt.gca().set_xlabel("frequency", y=1.1)
plt.gca().set_xlim([0,10])

(0.0, 10.0)